/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file
    Definitions of the sfrhist_snapshot methods.
    \ingroup sfrhist */

// $Id$

#include "sfrhist_snapshot.h"
#include <set>
#include "boost/lambda/lambda.hpp"
#include "boost/lambda/bind.hpp"
#include "boost/function.hpp"
#include "boost/lexical_cast.hpp"
#include "boost/foreach.hpp"

#include "mcrx-types.h"
#include "mcrx-debug.h"
#include "blitz-fits.h"
#include "model.h"
#include "fits-utilities.h"
#include "CCfits/CCfits"
#include "vecops.h"

using namespace std;
using namespace CCfits;
using namespace blitz;

/** Constructor loads the snapshot as well as the progenitor
    snapshots. The preferences object should contain the info from the
    GADGET parameter file, specifically the unit and softening length
    keywords. */
mcrx::sfrhist_snapshot::sfrhist_snapshot (const Preferences& gp,
					  const Preferences& p,
					  const string& snapshot_file, 
					  const string& snapshot_name,
					  const vector<string>& progenitors) : 
  Snapshot (gp, p, snapshot_file, snapshot_name)
{
  using namespace boost::lambda;

  // immediately convert our own units
  makephysicalunits();

  // things that only make sense in the sfrhist context, using the
  // parameter file preferences object and not the gadget paramater
  // prefs.

  // apply the star_radius_factor keyword
  const double star_radius_factor = p.defined("star_radius_factor") ?
    p.getValue("star_radius_factor", double()) : 1.0;
  cout << "Multiplying star particle radius by " << star_radius_factor << endl;
  BOOST_FOREACH(species s, star_species() ) {
    transform(set(s).h_.begin(), set(s).h_.end(),
	      set(s).h_.begin(),
	      _1 * star_radius_factor);
  }
  // do same for bh_radius_factor keyword
  const double bh_radius_factor = p.defined("bh_radius_factor") ?
    p.getValue("bh_radius_factor", double()) : 1.0;
  cout << "Multiplying BH particle radius by " << bh_radius_factor << endl;
  transform(bh_particles.h_.begin(), bh_particles.h_.end(),
	    bh_particles.h_.begin(),
	    _1 * bh_radius_factor);


  // Load progenitor snapshots
  for (int i = 0; i < progenitors.size(); ++i){
    boost::shared_ptr<Snapshot> is (new Snapshot (gp, p, progenitors[i],""));
    is->makephysicalunits(); 
    isnaps.push_back(is);
  }
  // snapshot class should make sure particle numbers are consistent with the progenitors

  min_id=1;
  if(p.getValue("nbodycod",string())=="arepo") {
    min_id=blitz::huge(int());
    for(int s=0; s < n_species; ++s)
      if(sets()[s]->n()>0)
	min_id=std::min(min_id, sets()[s]->min_ID());
    cout << "Reading Arepo snapshot, min ID is " << min_id << endl;
  }

  // Set up vector of progenitor particle number
  transform(isnaps.begin(), isnaps.end(), back_inserter(isnap_sizes), 
	    bind(&Snapshot::ntot,boost::cref(*_1)));

  // Allocate SED vectors (not the arrays themselves though)
  SED_set.resize(n_species);
  L_bol_set.resize(n_species);
}


/** Calculates the creation masses of the old populations based on mass
    loss info from the stellar model. */
void mcrx::sfrhist_snapshot::calculate_oldpop_massloss(const stellarmodel& m)
{
  disk_particles.calculate_creation_mass(time(), m);
  bulge_particles.calculate_creation_mass(time(), m);
}


/** Assigns the SED of the stellar particles from the model. Because
 the SED will potentially depend not just on the properties of the
 individual particles but on their aggregate properties (with
 clustering of MAPPINGS particles, for example), the particles are put
 in a list and the entire list sent to the model for SED assignment.
 The model has access to the full particle data and should return the
 physical SED of the particle, INCLUDING mass normalization. */
void
mcrx::sfrhist_snapshot::calculate_stellar_SEDs(const stellarmodel& m)
{
  // NOTE: there's the mass loss issue with old stars. This must be
  // handled before this point. It's really a matter of setting the
  // m_creation field of those species correctly and is not even
  // really the responsibility of this class...  The problem is that
  // the info for doing this comes from the stellar model, but there's
  // not really a problem making a special function for this in the
  // Snapshot class itself.

  // check that units are consistent
  check_consistency(units(), m.units(), true);
  
  cout << "Calculating stellar SEDs" << endl;

  vector<star_particle_ptr> particles;
 
  // Loop over stellar species
  for (std::set<species>::const_iterator s=star_species().begin();
       s != star_species().end(); s++) {

    // syntactic conveniences
    star_particle_set& current_set = 
      *dynamic_cast<star_particle_set*>(sets()[*s]);
    const int npart = current_set.n();
    
    // Allocate SED
    SED_set[*s].resize(npart, m.n_lambda());
    
    // push the particles onto the vector
    for (int i=0; i < npart; ++i) {
      particles.push_back(star_particle_ptr(current_set, i, SED_set[*s]));
    }
  }

  // and ask the model to calculate SEDs. If the age is NaN then the
  // particle should be black. Note that the model should set the
  // actual SED with correct mass normalization
  m.calculate_SEDs(particles, time());

  // we need the following info from the model for later when we save the file
  // we need this to know how wide to make the table
  n_lambda = m.n_lambda(); 
  SED_unit =m.units().get("luminosity")+"/"+ m.units().get("wavelength");
  L_unit = m.units().get("luminosity");
}



/** Assigns the SED of the black hole particles from the bh model. This
    function simply asks the bh model for the SED based on the info in
    the particle. The model has access to the full particle data and
    should return the physical SED of the particle. NOTE: This must be
    run AFTER calculate_stellar_SEDs has been run and the BH SEDs have
    been resampled onto the wavelengths used for the stellar calculation.*/
void
mcrx::sfrhist_snapshot::calculate_BH_SEDs(const bhmodel& m)
{
  // check that units are consistent
  check_consistency (units (), m.units(), true);
  
  // Check that the number of wavelengths is correct.
  if (n_lambda != m.n_lambda() ) {
     cerr << "FATAL ERROR: Wavelengths are not consistent between stellar model and black hole model!" << endl;
     exit (1);
  }

  cout << "Calculating BH SEDs" << endl;

  const int npart = bh_particles.n();

  // Allocate SED.
  SED_set[bh].resize(npart, m.n_lambda());

  // Loop over particles
  for (int i=0; i < npart; ++i) {
    // and get SED.
    SED_set[bh](i, blitz::Range::all ()) = m.get_SED (bh_particles, i);
  }

}




/** Deduces which object and ID number in the initial conditions the
    particle had. For new stellar particles the ID of the parent gas
    particle is returned. This is strictly a gadget-specific function
    so at some point we might need to make this mechanism more
    generic. If there are no isnaps defined, it returns (-1, ID).
    Note that this function relies on the fact that in the initial
    snapshots, ID numbers are contiguous, so not loading any of the
    particle sets using the "load_xxx false" keyword will make this
    function fail. Note that this function does not work for Arepo
    mergers because mesh cells do not retain their identity like
    particles do. */
pair<int, unsigned int>
mcrx::sfrhist_snapshot::deduce_progenitor(unsigned int ID) const
{
  if (isnap_sizes.size()==0) {
    // no isnaps defined
    return make_pair(-1, ID);
  }

  unsigned int parent_snap_first = 0;
  int parent_snap = 0;
  assert(min_id==1 || isnap_sizes.size()<=1);
  while (parent_snap < isnap_sizes.size()) {
    // see if ID came from this snapshot
    const int ntot = isnap_sizes [parent_snap];
    if ((ID - min_id - parent_snap_first)< ntot) // ID numbering starts w/min_id hence <
      // success!
      return make_pair (parent_snap, ID - parent_snap_first);

    // not in this snapshot, increment and try next one
    parent_snap_first += ntot;
    ++parent_snap;
  }

  // if we come out here, the particle ID is greater than the sum of
  // the number of particles in all initial conditions, which means it
  // is a new star particle or a gas particle that has spawned a new star.
  // We then pick up the parent ID and retry on that. (ID_index will
  // throw if we've screwed up.)

  // Try first to see if it is a gas particle
  try {
    const unsigned int newID = gas_particles.pid_[gas_particles.find_ID(ID)];
    // make sure we don't go into an infinite loop
    assert(newID != ID);
    return deduce_progenitor (newID);
  }
  catch (particle_set_generic::illegal_id&) {};

  // It could also be a second-generation new star particle
  try {
    const unsigned int newID = nstar_particles.pid_[nstar_particles.find_ID(ID)];
    assert(newID != ID);
    return deduce_progenitor (newID);
  }
  catch (particle_set_generic::illegal_id&) {};

  // If we haven't found it now it's a bad ID, so we throw.
  throw particle_set_generic::illegal_id();

}


/** Assigns age (actually formation time) for the particles of the
    specified species from the functors in the vector. The entries in
    the vector are for the different progenitors. The function
    determines the progenitor for each particle, calls that functor
    and sets the metallicity and age for the particle to the result.
    The functors get 4 arguments: the index of the current particle, a
    pointer to the current particle set, the index of the progenitor
    particle, and a pointer to the progenitor particle set. The
    functor should return a pair (formation time, metallicity) of the
    particle. Typically, the functor will want to ADD the current
    metallicity to the one determined from the progenitor position,
    this function does NOT do that. For new star particles, the
    formation time is already known and the functor should just return
    that value if it wants to use it. (Again, it can set the formation
    time as it wants.)  Returning NaN for age is a special case which
    means the particle has no emission, this can be used for "absent"
    populations. For gas species, the tform value is ignored.  */
void
mcrx::sfrhist_snapshot::assign_age_metallicity (const species s,
						const vector<boost::shared_ptr<T_age_metallicity_functor> >& f_age_metallicity)
{
  // Because neither age nor metallicity are part of the generic
  // particle sets, we need to extract those vectors. We just switch
  // on the types and dynamic_cast the species to the selected
  // types. This is not very extensible but should not be a big deal
  // here. A bad_cast will be thrown if you try to set the metallicity
  // of e.g. the DM particles, but that doesn't seem like a bad
  // thing...
  vector<T_float>* species_tf = 0;
  vector<T_float>* species_z = 0;

  if (s == gas) {
    // species is gas, set metallicity vector pointer and keep tf null.
    species_z = &(dynamic_cast<gas_particle_set*>(sets()[s])->z_);
  }
  else if (star_species().find(s) != star_species().end()) {
    // species is star, set both pointers.
    species_tf = &(dynamic_cast<star_particle_set*>(sets()[s])->tf_);
    species_z = &(dynamic_cast<star_particle_set*>(sets()[s])->z_);
  }

  // Loop over particles 
  for (int i=0; i< n(s); ++i) {
    
    // get the "progenitor" ID and object
    const int this_id = sets()[s]->id(i);
    const pair<int, unsigned int> p = deduce_progenitor (this_id);
    const int prog_obj = p.first;

    pair<T_float, T_float> age_z;

    if (prog_obj<0) {
      // this means no progenitor was found, so the particle is an
      // orphan that does not belong to any galaxy. In that case we
      // use the particle itself. (It only makes sense to define the
      // orphaned particles in a way that is spatially independent,
      // since we have no spatial info.)
      age_z = 
	(*(f_age_metallicity[nobject()]))(i, *sets()[s],
					  i, *sets()[s]);
    }
    else {
      const unsigned int prog_id = p.second;

    // there is a subtlety here, we now know which object the
    // progenitor is in, and what its ID is. However, we don't know
    // what species it's in, because we might have a new star particle
    // in which case the progenitor is a gas particle. This means we
    // must use the most generic search function that looks through
    // all the species in the progenitor snapshot.

    // find particle vector index in progenitor snapshot
    const pair<species, int> prog_species_index = 
      isnaps[prog_obj]->find_ID(prog_id);
    const species prog_species = prog_species_index.first;
    const int prog_index = prog_species_index.second;

    // call functor to get age and metallicity.
    age_z = 
      (*(f_age_metallicity[prog_obj]))(i, *sets()[s], prog_index, 
				       *isnaps[prog_obj]->sets()[prog_species]);
    }
    // and assign (check for null for tform)
    // (actually, the metallicity is ADDED to the metallicity from the snapshot)
    (*species_z)[i] += age_z.second;
    if (species_tf)
      (*species_tf)[i] = age_z.first;
  }
}


/** Adds the center of mass (actually integrates both m_i*pos_i and
    m_i) for the particles of species s, keeping track of objects
    separately and subject to an inclusion condition. Note that cm and
    m are changed. */
void 
mcrx::sfrhist_snapshot::
add_cm(species s,
       boost::function<bool(int,const particle_set_generic*const,int)> predicate,
       vector<vec3d>& cm, vector<double>& m) const
{
  for(int i=0; i< n(s); ++i) {
    pair<int, unsigned int> p = deduce_progenitor(sets()[s]->id(i));
    int prog_object = p.first;
    if(prog_object<0)
      // particle does not belong to an object. put it in the generic
      // category which is the last entry
      prog_object=cm.size()-1;

    if(predicate(prog_object, sets()[s], i)) {
      const T_float current_m = sets()[s]->m(i);
      m[prog_object]+= current_m;
      cm[prog_object]+= sets()[s]->pos(i)*current_m;
    }
  }
}

/** There is some resolution problem with the - operator and lambda,
    so defining this is necessary. */
mcrx::vec3d vec3d_minus(const mcrx::vec3d& a, const mcrx::vec3d& b) {return a-b;}

/** Calculates the centers of the component galaxies by calculating
    the cm in apertures of decreasing size. */
vector<mcrx::vec3d> mcrx::sfrhist_snapshot::calculate_centers () const
{
  cout << "Calculating centers of galaxies:" << endl;

  using namespace boost::lambda;
  using namespace boost;

  double aperture=1e6;
  double end_aperture=aperture;
  if(n(disk))
    end_aperture =min(end_aperture,
		      *min_element(disk_particles.h_.begin(),
				   disk_particles.h_.end()));
  if(n(bulge))
    end_aperture =min(end_aperture,
		      *min_element(bulge_particles.h_.begin(),
				   bulge_particles.h_.end()));
  if(n(nstar))
    end_aperture =min(end_aperture,
		      *min_element(nstar_particles.h_.begin(),
				   nstar_particles.h_.end()));

  // make vector of particle number of initial snapshots
  const int n=isnaps.size()+1;  
  vector<int> max_particle(n);
  vector<vec3d> cm(n);

  bool not_done=true;
  bool first=true;
  T_float convergence_factor=1;

  typedef const vec3d&(particle_set_generic::*T_posfn)(int)const;
  T_posfn posfunc=&particle_set_generic::pos;
  // the following as a lambda expr that says: 
  // mag(_2.pos(_3) - cm[_1] < aperture,
  // ie true if particle is within aperture of the current cm
  const boost::function<bool(int,const particle_set_generic*const,int)> predicate =
    bind(&mag, 
	 bind(vec3d_minus,
	      bind(posfunc,_2,_3),
	      var(cm)[_1]))
    < var(aperture);
    
  while(true) {
    vector<vec3d> new_cm(n, vec3d(0,0,0));
    vector<double> m(n, 0.);
    add_cm(disk, predicate, new_cm, m);
    add_cm(bulge, predicate, new_cm, m);
    add_cm(nstar, predicate, new_cm, m);
    bool valid=true;
    for(int i=0; i<n; ++i) {
      new_cm[i]=new_cm[i]/m[i];
      // we don't include the "no galaxy" particles in validity unless 
      // all particles are "no galaxy"
      valid&=  (all(new_cm[i]==new_cm[i]) && (m[i]>0) &&
		all(abs(new_cm[i])<blitz::huge(T_float()))) ||
	((n>1) && (i==n-1) && (m[i]==0));
      DEBUG(1,cout << " Galaxy " << i << " center: " << new_cm[i] << '\t' << m[i] << '\t' << valid << endl;);
    }
    DEBUG(1,cout << "Aperture " << aperture << '\t' << convergence_factor << endl;);
    if (!valid) {
      // we did not converge, maybe no particles in aperture? make
      // convergence factor smaller and try again
      aperture*=(1+convergence_factor);
      convergence_factor/=2;
      aperture/=(1+convergence_factor);
      if(convergence_factor<1e-3) {
	// something's wrong, abort
	cerr << "Could not converge to a good solution, aborting\n";
	break;
      }
      DEBUG(1,cout << "Step failed, backing up\n";);
      continue;
    }

    cm=new_cm;
    aperture/=(1+convergence_factor);
    convergence_factor*=1.25;
    // iterate until aperture is smaller than the softening of the particles
    if(aperture<end_aperture)
      break;
  }

  // when we come out here, we have converged centers (hopefully) in cm
  // return them

  for(int i=0; i<n-1; ++i)
    cout << " Galaxy " << i << " center: " << cm[i] << endl;
  if(all(cm.back()==cm.back())) {
    // only print if there are such particles, if there aren't cm is NaN.
    cout << "Particles not belonging to a galaxy: " << cm[n-1] << endl;
  }
  else {
    // there are no unbound paticles, delete that entry.
    cm.pop_back();
  }

  return cm;
}
  

/** Changes the origin of the particles. */
void mcrx::sfrhist_snapshot::change_origin(const vec3d& delta)
{
  using namespace boost::lambda;
  for (int s=0; s<n_species; ++s)
    transform(sets()[s]->pos_.begin(), sets()[s]->pos_.end(), 
	      sets()[s]->pos_.begin(),
	      _1-=constant(delta));
}


/** Calculates the bolometric luminosity of the particles by
    integrating the SEDs over wavelength. */
void mcrx::sfrhist_snapshot::calculate_L_bol (const stellarmodel& m) 
{
  using namespace blitz;
  cout << "Calculating bolometric luminosities" << endl;
  
  // This is done generically for all species
  for (int s=0; s<n_species; ++s) {
    L_bol_set[s].resize(n(species(s)));
    if (SED_set[s].size() > 0)
      for(int i=0; i<SED_set[s].extent(firstDim); ++i) {
	L_bol_set[s](i) = integrate_quantity(SED_set[s](i, Range::all()),
					     m.lambda(),
					     true);
	assert(L_bol_set[s](i)<1e60);
      }
    else
      L_bol_set[s] = 0.;
  }
}

void
mcrx::sfrhist_snapshot::write_fits_table (const string& output_file_name) const
{
  const string hdu_name = "PARTICLEDATA";

  // Check if we have physical units
  const bool physical_units =(massunit_ == 1) && (lengthunit_ == 1) &&
    (timeunit_ == 1);
  if (!physical_units) {
    cerr << "Error: Snapshot is not in default units." << endl;
    exit(1);
  }

  // Create FITS table.
  auto_ptr<FITS> output_file;
  output_file.reset(new FITS (output_file_name, Write ));

  {
    Table* output = output_file->addTable(hdu_name,0);
    output->addColumn(Tuint, "ID", 1, "");
    output->addColumn(Tdouble, "mass", 1, units().get("mass") );
    output->addColumn(Tdouble, "position", 3, units().get("length"));
    output->addColumn(Tdouble, "velocity", 3, 
		      units().get("length")+"/"+units().get("time"));
    output->addColumn(Tdouble, "radius", 1,units().get("length"));
    output->addColumn(Tdouble, "metallicity", 1, "");
    output->addColumn(Tdouble, "formation_time", 1, units().get("time"));
    output->addColumn(Tint, "parent_ID", 1, "");
    output->addColumn(Tdouble, "creation_mass", 1, units().get("mass"));
    
    output->addColumn(Tdouble, "age", 1, units().get("time"));
    output->addColumn(Tfloat , "L_lambda", n_lambda, SED_unit);
    output->addColumn(Tdouble, "L_bol", 1, L_unit);
    
    output->addKey("logflux", true,
		   "Column L_lambda values are log (L_lambda)");
  }

  // Loop over stellar species
  array_1 L_lambda_tot(n_lambda);
  L_lambda_tot = 0;
  for (std::set<species>::const_iterator s=star_species().begin();
       s != star_species().end(); s++) {

    // Write particles. The reason for repeated reopening of the file
    // is that writing is done with cfitsio, so CCfits will get
    // confused if we don't reopen the file every time.
    output_file.reset(new FITS (output_file_name, Write ));
    write_fits_set(open_HDU(*output_file, hdu_name), *s);

    L_lambda_tot += sum(SED_set[*s](tensor::j, tensor::i), tensor::j);
  }

  // write BH partcles
  output_file.reset(new FITS (output_file_name, Write ));
  write_fits_set(open_HDU(*output_file, hdu_name), bh);
  L_lambda_tot += sum(SED_set[bh](tensor::j, tensor::i), tensor::j);

  // also write INTEGRATED_QUANTITIES HDU
  output_file.reset(new FITS (output_file_name, Write ));
  const string iq_hdu_name = "INTEGRATED_QUANTITIES";
  cout << "Writing integrated grid SED to HDU " << iq_hdu_name << endl;
  ExtHDU* iq_hdu;
  try {
      iq_hdu = &open_HDU(*output_file,iq_hdu_name);
    }
  catch(...) {
    iq_hdu = output_file->addTable(string (iq_hdu_name) , 0);
  }
  iq_hdu->writeComment("This HDU contains the integrated SED of all particles (even those outside of the grid)");
  iq_hdu->addColumn(Tdouble, "L_lambda", 1, SED_unit );
  write (iq_hdu->column("L_lambda" ), L_lambda_tot, 1 );
}


/** Writes the data for a (star) particle set to a FITS table.  The only
    complicated thing here is what we are going to do with the BH
    particles... */
void
mcrx::sfrhist_snapshot::write_fits_set (CCfits::ExtHDU& output, 
					species s) const
{
  using namespace CCfits;
  using namespace boost::lambda;

  CCfits::Column& c_age = output.column("age") ;
  CCfits::Column& c_L_bol = output.column("L_bol") ;
  CCfits::Column& c_L_lambda = output.column("L_lambda") ;

  // continue filling table at already existing row count
  const int start_row = c_age.rows()+1;
  const particle_set_generic& this_set = *sets()[s];
  const int n_rows = this_set.n();
  
  // We don't write formation time but rather age, so calculate it
  // here (if the particles have ages, otherwise write NaN)
  vector<T_float> age;
  try {
    const star_particle_set& this_set = 
      dynamic_cast<const star_particle_set&>(*sets()[s]);
    transform(this_set.tf_.begin(), this_set.tf_.end(), back_inserter(age),
	      time() - _1);
    // for Gadget hdf5 files, the snaptime is saved as a double but
    // tform as a float. This means that particles that form
    // immediately before snapshot writing can end up with tforms
    // larger than snaptime. in this case we just set the age to zero
    replace_if(age.begin(), age.end(), _1 < 0.0, 0.0);
  }
  catch (std::bad_cast) {
    age.assign(n_rows, blitz::quiet_NaN(T_float()));
  }

  // check that SED is defined
  array_2 current_sed;
  if(SED_set[s].size()==0) {
    cout << "Warning: No SED defined for species " << s
	 << ", assuming zero" << endl;
    current_sed.resize(n_rows, c_L_lambda.repeat() );
    current_sed=0;
  }
  else {
    current_sed.reference(SED_set[s]);
    ASSERT_ALL(current_sed<1e60);
    ASSERT_ALL(L_bol_set[s]<1e60);
  }

  // to get reasonable I/O performance from cfitsio we need to write
  // in chunks.  
  const long n_chunk =  output.getRowsize();
  cout << "      chunk size is " << n_chunk << " rows" << endl;
  int j = 0;
  int status=0;
  while (j < n_rows ) {
    const int this_chunk = (n_rows- j > n_chunk)?  
      n_chunk: n_rows-j;
    const int current_row = j+ start_row;

    // The star_particle fields are written by the set itself
    this_set.write_table_chunk(current_row, j, this_chunk, output);

    // write age
    make_onedva(age).write_column(output, c_age, current_row, 
				  j, this_chunk);

    // Write the luminosity data
    make_onedva(L_bol_set[s]).write_column(output, c_L_bol, current_row,
					   j, this_chunk);

    // We Write log of the SED.  0 values are problematic so add a
    // very small number to all of them. 
    assert (j+ this_chunk - 1 <= current_sed.rows());
    const array_2 logSED (log10 (current_sed(Range (j ,j + this_chunk-1 ), 
					    Range::all ())+
				 blitz::tiny (T_float ())));
    write (output.fitsPointer(), c_L_lambda.index(), logSED, current_row );

    j+= this_chunk;
  }
  fits_report_error(stderr,status);
}
